* note -------------------------------------------------------------------------

* this is one of the two replication scripts for:
*
* umit, r., & schaffer, l. m. (2020). attitudes towards carbon taxes across 
* europe: the role of perceived uncertainty and self-interest. energy policy,
* 140, 1-7.
*
* the other replication script is called 02_analysis_r.R.
*
* running this script requires two third-party packages, which can be 
* installed by removing the leading asterics from the lines 22 and 23 below.
*
* resul umit, 19 july 2022

* declare version --------------------------------------------------------------
version 15.1


* install additional packages --------------------------------------------------

* ssc install estout
* ssc install mlt


* prepare the data -------------------------------------------------------------

* import the original ess dataset
use ESS8e02_1.dta, clear
 
* merge with context data
merge m:1 cntry using 03_data_context.dta
drop _merge

* adjust the income variable
tostring hinctnta, generate(decile_number)
gen merge_key = cntry + "_" + decile_number
merge m:1 merge_key using 04_data_income.dta
drop _merge

* recode the dependent variable so that higher values indicate an increase
recode inctxff 1=5 2=4 4=2 5=1

* label values of the dependent variable
label define support 1 "Strongly against" 2 "Somewhat against"                  ///
3 "Neither in favour nor against" 4 "Somewhat in favour" 5 "Strongly in favour"
label values inctxff support

* recode consumption as dependence
recode cflsenr 0=10 1=9 2=8 3=7 4=6 5=5 6=4 7=3 8=2 9=1 10=0

* label values of the the dependence variable
label define dependence 0 "Completely confident" 1 "Not at all confident"
label values cflsenr dependence

* recode other as missing for education
replace eisced =. if eisced > 7

* recode gender so that male=0 and female=1
recode gndr 1=0 2=1

* label the values of the gender variable
label define gender 0 "Male" 1 "Female"
label values gndr gender

* recode gdp in 10,000 dollars
replace gdp = gdp / 10000

* recode "cannot reduce" as missing for rdcenr
replace rdcenr =. if rdcenr == 55   

* create the political trust variable
egen trstpol=rmean(trstprl trstplt trstprt)

* keep the variables needed for analysis
keep inctxff trstpol psppsgva cflsenr domicil                                   ///
agea eisced gndr hinctnta adj_income lrscale wrclmch wrenexp eneffap rdcenr     ///
carbon_tax party_position gdp cntry pspwght pweight

* label the variables
label variable inctxff        "Support for Taxes"

label variable trstpol        "Political Trust"
label variable psppsgva       "Political Efficacy"
label variable cflsenr        "Energy Dependence"
label variable domicil        "Rural Area"

label variable agea           "Age"
label variable eisced         "Education"
label variable gndr           "Female"
label variable adj_income     "Household Income"
label variable lrscale        "Left-Right Orientation"
label variable wrclmch        "Climate Worries"
label variable wrenexp        "Cost Worries"
label variable eneffap        "Efficiency"
label variable rdcenr         "Curtailment"

label variable carbon_tax     "Carbon Tax"
label variable party_position "Party Position"
label variable gdp            "GDP"

label variable cntry          "Country"

* save the clean dataset
save data_ess.dta, replace


* table 1 ----------------------------------------------------------------------

* estimate the models
xtmixed inctxff trstpol psppsgva cflsenr domicil || cntry:
estimates store model1
mltrsq

xtmixed inctxff trstpol psppsgva cflsenr domicil                                ///
agea eisced gndr adj_income lrscale wrclmch wrenexp eneffap rdcenr              ///
|| cntry:
estimates store model2
mltrsq

xtmixed inctxff trstpol psppsgva cflsenr domicil                                ///
carbon_tax party_position gdp || cntry:
estimates store model3
mltrsq

xtmixed inctxff trstpol psppsgva cflsenr domicil                                ///
agea eisced gndr adj_income lrscale wrclmch wrenexp eneffap rdcenr              ///
carbon_tax party_position gdp || cntry:
estimates store model4
mltrsq

* print the table
estout model1 model2 model3 model4,                                             ///
cells(b(star fmt(3)) se(par fmt(3))) label nonumber eqlabel("")                 ///
keep(trstpol psppsgva cflsenr domicil cflsenr  _cons)                           ///
varlabels(_cons Constant)                                                       ///
title("Multilevel linear regression models---summary results")


* estimate predicted probabilities for figure 3 --------------------------------

* political trust
** re-estimate the model
quietly xtmixed inctxff c.trstpol psppsgva cflsenr domicil                      ///
agea eisced gndr adj_income lrscale wrclmch wrenexp eneffap rdcenr              ///
carbon_tax party_position gdp || cntry:
** estimate the predictions
estpost margins, at(trstpol=(0(1)10))
** save the results, to be plotted in r
esttab using predictions_linear_trust.csv, cell("b ci_l ci_u") plain replace

* political efficacy
** re-estimate the model
quietly xtmixed inctxff c.trstpol psppsgva cflsenr domicil                      ///
agea eisced gndr adj_income lrscale wrclmch wrenexp eneffap rdcenr              ///
carbon_tax party_position gdp || cntry:
** estimate the predictions
estpost margins, at(psppsgva=(1(1)5))
** save the results, to be plotted in r
esttab using predictions_linear_efficacy.csv, cell("b ci_l ci_u") plain replace

* energy dependence
** re-estimate the model
quietly xtmixed inctxff c.trstpol psppsgva cflsenr domicil                      ///
agea eisced gndr adj_income lrscale wrclmch wrenexp eneffap rdcenr              ///
carbon_tax party_position gdp || cntry:
** estimate the predictions
estpost margins, at(cflsenr=(0(1)10))
** save the results, to be plotted in r
esttab using predictions_linear_dependence.csv, cell("b ci_l ci_u") plain       ///
replace

* rural areas
** re-estimate the model
quietly xtmixed inctxff c.trstpol psppsgva cflsenr domicil                      ///
agea eisced gndr adj_income lrscale wrclmch wrenexp eneffap rdcenr              ///
carbon_tax party_position gdp || cntry:
** estimate the predictions
estpost margins, at(domicil=(1(1)5))
** save the results, to be plotted in r
esttab using predictions_linear_rural.csv, cell("b ci_l ci_u") plain replace

* table s1 ---------------------------------------------------------------------

eststo: quietly estpost tabstat inctxff trstpol psppsgva cflsenr domicil        ///
agea eisced gndr adj_income lrscale wrclmch wrenexp eneffap rdcenr              ///
carbon_tax party_position gdp,                                                  ///
statistics(count mean median sd min max) columns(statistics)

esttab,                                                                         ///
cells("count mean(fmt(2)) p50(fmt(2)) sd(fmt(2)) min(fmt(2)) max(fmt(2))")      ///
noobs label nonumber title("Descriptive statistics")
eststo clear

* table s3 ---------------------------------------------------------------------

estout model1 model2 model3 model4,                                             ///
cells(b(star fmt(3)) se(par fmt(3))) label nonumber eqlabel("")                 ///
varlabels(_cons Constant)                                                       ///
title("Multilevel linear regression models---complete results")


* table s4 ---------------------------------------------------------------------

* estimate the models
xtmixed inctxff trstpol psppsgva cflsenr domicil || cntry:
estimates store modela1
mltrsq

xtmixed inctxff trstpol psppsgva cflsenr domicil                                ///
agea || cntry:
estimates store modela2
mltrsq

xtmixed inctxff trstpol psppsgva cflsenr domicil                                ///
agea eisced || cntry:
estimates store modela3
mltrsq

xtmixed inctxff trstpol psppsgva cflsenr domicil                                ///
agea eisced gndr || cntry:
estimates store modela4
mltrsq

xtmixed inctxff trstpol psppsgva cflsenr domicil                                ///
agea eisced gndr adj_income || cntry:
estimates store modela5
mltrsq

xtmixed inctxff trstpol psppsgva cflsenr domicil                                ///
agea eisced gndr adj_income lrscale || cntry:
estimates store modela6
mltrsq

xtmixed inctxff trstpol psppsgva cflsenr domicil                                ///
agea eisced gndr adj_income lrscale wrclmch || cntry:
estimates store modela7
mltrsq

xtmixed inctxff trstpol psppsgva cflsenr domicil                                ///
agea eisced gndr adj_income lrscale wrclmch wrenexp || cntry:
estimates store modela8
mltrsq

xtmixed inctxff trstpol psppsgva cflsenr domicil                                ///
agea eisced gndr adj_income lrscale wrclmch wrenexp eneffap || cntry:
estimates store modela9
mltrsq

xtmixed inctxff trstpol psppsgva cflsenr domicil                                ///
agea eisced gndr adj_income lrscale wrclmch wrenexp eneffap rdcenr              ///
|| cntry:
estimates store modela10
mltrsq

xtmixed inctxff trstpol psppsgva cflsenr domicil                                ///
agea eisced gndr adj_income lrscale wrclmch wrenexp eneffap rdcenr              ///
carbon_tax || cntry:
estimates store modela11
mltrsq

xtmixed inctxff trstpol psppsgva cflsenr domicil                                ///
agea eisced gndr adj_income lrscale wrclmch wrenexp eneffap rdcenr              ///
carbon_tax party_position || cntry:
estimates store modela12
mltrsq

xtmixed inctxff trstpol psppsgva cflsenr domicil                                ///
agea eisced gndr adj_income lrscale wrclmch wrenexp eneffap rdcenr              ///
carbon_tax party_position gdp || cntry:
estimates store modela13
mltrsq

* print the table
estout modela1 modela2 modela3 modela4 modela5 modela6 modela7 modela8 modela9  ///
modela10 modela11 modela12 modela13 ,                                           ///
cells(b(star fmt(3)) se(par fmt(3))) label nonumber eqlabel("")                 ///
varlabels(_cons Constant)                                                       ///
title("Multilevel linear regression models---adding one control at a time")


* estimate vif scores for table s5 ---------------------------------------------
quietly regress inctxff trstpol psppsgva cflsenr domicil                        ///
agea eisced gndr adj_income lrscale wrclmch wrenexp eneffap rdcenr              ///
carbon_tax party_position gdp

estat vif


* table s6 ---------------------------------------------------------------------

* estimate the models
meologit inctxff trstpol psppsgva cflsenr domicil || cntry:
estimates store model5

meologit inctxff trstpol psppsgva cflsenr domicil                               ///
agea eisced gndr adj_income lrscale wrclmch wrenexp eneffap rdcenr              ///
|| cntry:
estimates store model6

meologit inctxff trstpol psppsgva cflsenr domicil                               ///
carbon_tax party_position gdp || cntry:
estimates store model7

meologit inctxff c.trstpol psppsgva cflsenr domicil                             ///
agea eisced gndr adj_income lrscale wrclmch wrenexp eneffap rdcenr              ///
carbon_tax party_position gdp || cntry:
estimates store model8

* print the table
estout model5 model6 model7 model8,                                             ///
cells("b(star fmt(3)) se(par fmt(3))") label nonumber eqlabel("")               ///
title("Ordered logistic regression models")


* estimate predicted probabilities for figure s1 -----------------------------

* political trust
** re-estimate the model
quietly meologit inctxff c.trstpol psppsgva cflsenr domicil                     ///
agea eisced gndr adj_income lrscale wrclmch wrenexp eneffap rdcenr              ///
carbon_tax party_position gdp || cntry:
** estimate the predictions
estpost margins, at(trstpol=(0(1)10)) expression(predict(pr outcome(4))         ///
+ predict(pr outcome(5)))
** save the results, to be plotted in r
esttab using predictions_ordinal_trust.csv, cell("b ci_l ci_u") plain replace

* political efficacy
** re-estimate the model
quietly meologit inctxff c.trstpol psppsgva cflsenr domicil                     ///
agea eisced gndr adj_income lrscale wrclmch wrenexp eneffap rdcenr              ///
carbon_tax party_position gdp || cntry:
** estimate the predictions
estpost margins, at(psppsgva=(1(1)5)) expression(predict(pr outcome(4))         ///
+ predict(pr outcome(5)))
** save the results, to be plotted in r
esttab using predictions_ordinal_efficacy.csv, cell("b ci_l ci_u") plain replace

* energy dependence
** re-estimate the model
quietly meologit inctxff c.trstpol psppsgva cflsenr domicil                     ///
agea eisced gndr adj_income lrscale wrclmch wrenexp eneffap rdcenr              ///
carbon_tax party_position gdp || cntry:
** estimate the predictions
estpost margins, at(cflsenr=(0(1)10)) expression(predict(pr outcome(4))         ///
+ predict(pr outcome(5)))
** save the results, to be plotted in r
esttab using predictions_ordinal_dependence.csv, cell("b ci_l ci_u") plain      ///
replace

* rural areas
** re-estimate the model
quietly meologit inctxff c.trstpol psppsgva cflsenr domicil                     ///
agea eisced gndr adj_income lrscale wrclmch wrenexp eneffap rdcenr              ///
carbon_tax party_position gdp || cntry:
** estimate the predictions
estpost margins, at(domicil=(1(1)5))  expression(predict(pr outcome(4))         ///
+ predict(pr outcome(5)))
** save the results, to be plotted in r
esttab using predictions_ordinal_rural.csv, cell("b ci_l ci_u") plain replace
